home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 2 / AACD 2.iso / AACD / Sound / LAME / Source / l3psy.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-06-03  |  36.8 KB  |  1,409 lines
  1. /**********************************************************************
  2.  * ISO MPEG Audio Subgroup Software Simulation Group (1996)
  3.  * ISO 13818-3 MPEG-2 Audio Encoder - Lower Sampling Frequency Extension
  4.  *
  5.  * $Id: l3psy.c,v 1.2 1998/10/05 17:06:48 larsi Exp $
  6.  *
  7.  * $Log: l3psy.c,v $
  8.  * Revision 1.2  1998/10/05 17:06:48  larsi
  9.  * *** empty log message ***
  10.  *
  11.  * Revision 1.1.1.1  1998/10/05 14:47:18  larsi
  12.  *
  13.  * Revision 1.2  1997/01/19 22:28:29  rowlands
  14.  * Layer 3 bug fixes from Seymour Shlien
  15.  *
  16.  * Revision 1.1  1996/02/14 04:04:23  rowlands
  17.  * Initial revision
  18.  *
  19.  * Received from Mike Coleman
  20.  **********************************************************************/
  21. /**********************************************************************
  22.  *   date   programmers         comment                               *
  23.  * 2/25/91  Davis Pan           start of version 1.0 records          *
  24.  * 5/10/91  W. Joseph Carter    Ported to Macintosh and Unix.         *
  25.  * 7/10/91  Earle Jennings      Ported to MsDos.                      *
  26.  *                              replace of floats with FLOAT          *
  27.  * 2/11/92  W. Joseph Carter    Fixed mem_alloc() arg for "absthr".   *
  28.  * 3/16/92  Masahiro Iwadare    Modification for Layer III            *
  29.  * 17/4/93  Masahiro Iwadare    Updated for IS Modification           *
  30.  **********************************************************************/
  31.  
  32. #include "common.h"
  33. #include "globalflags.h"
  34. #include "encoder.h"
  35. #include "l3psy.h"
  36. #include "l3side.h"
  37. #include <assert.h>
  38. #include "gtkanal.h"
  39. #define MSFREQ 20
  40.  
  41. #define maximum(x,y) ( (x>y) ? x : y )
  42. #define minimum(x,y) ( (x<y) ? x : y )
  43.  
  44. /* some different types of adaptive window switching */
  45.  
  46.  
  47. /* This mode will turn on short_blocks if there is a localized surge in energy */
  48. #define ENER_AWS  
  49.  
  50. /* from Gabriel Bouvigne <bouvigne@multimania.com> */
  51. #define AWS    
  52.  
  53. static int switch_pe=1800;
  54.  
  55.  
  56. void sprdngf1(FLOAT *, double *);
  57. void sprdngf2(double *, FLOAT *);
  58. static double s3_l[CBANDS][CBANDS]; /* needed global static by sprdngfs */
  59.  
  60. void L3para_read( double sfreq, int numlines[CBANDS], int partition_l[HBLKSIZE],
  61.           double minval[CBANDS], double qthr_l[CBANDS], double norm_l[CBANDS],
  62.           double s3_l[CBANDS][CBANDS], int partition_s[HBLKSIZE_s], double qthr_s[CBANDS_s],
  63.           double norm_s[CBANDS_s], double SNR_s[CBANDS_s],
  64.           int cbw_l[SBMAX_l], int bu_l[SBMAX_l], int bo_l[SBMAX_l],
  65.           double w1_l[SBMAX_l], double w2_l[SBMAX_l],
  66.           int cbw_s[SBMAX_s], int bu_s[SBMAX_s], int bo_s[SBMAX_s],
  67.           double w1_s[SBMAX_s], double w2_s[SBMAX_s] );
  68.                                     
  69.  
  70.  
  71.  
  72.  
  73.  
  74.  
  75.  
  76.  
  77.  
  78.  
  79.  
  80.  
  81. void L3psycho_energy( short int *buffer, 
  82.     FLOAT energy[HBLKSIZE],
  83.     FLOAT ax[HBLKSIZE], FLOAT bx[HBLKSIZE],
  84.     FLOAT energy_s[3][HBLKSIZE_s],
  85.     FLOAT ax_s[3][HBLKSIZE_s], FLOAT bx_s[3][HBLKSIZE_s],
  86.      int chn,int gr , layer * info,
  87.      int check_ms_stereo, double *ms_ener_ratio)
  88. {
  89.   static short int savebuffer[2][1344];
  90.   static double ms_ener[2];
  91.   static int sync_flush,flush,syncsize;
  92.   static int init=0;
  93.   static FLOAT scalefac;
  94.   static FLOAT window_s[BLKSIZE_s];
  95.   static FLOAT window[BLKSIZE];
  96.  
  97. #ifdef HAVEGTK
  98.   static FLOAT energy_save[2][HBLKSIZE];
  99. #endif
  100.  
  101.   short int *savebuf;
  102.   int i,j,k,sblock;
  103.   
  104.   static FLOAT wsamp_r_int[2][BLKSIZE];
  105.   FLOAT wsamp_rs[256],*wsamp_r; 
  106.  
  107.  
  108.   wsamp_r=wsamp_r_int[chn]; 
  109.   savebuf = &savebuffer[chn][0];
  110.   
  111.   if(init==0){
  112.     init=1;
  113.     memset((char *) savebuffer, 0, sizeof(savebuffer));    
  114.     if (gpsycho) {sync_flush=WINDELAY; flush=576; syncsize=WINDELAY+576;}
  115.     else        {sync_flush=768; flush=576; syncsize=1344;} 
  116.     
  117.     scalefac=1.0;
  118.     if (force_ms) scalefac=sqrt(2.0);
  119.     
  120.     /* calculate HANN window coefficients */
  121.     /*   for(i=0;i<BLKSIZE;i++)  window[i]  =0.5*(1-cos(2.0*PI*i/(BLKSIZE-1.0)));*/
  122.     for(i=0;i<BLKSIZE;i++)  window[i]  =0.5*(1-cos(2.0*PI*(i-0.5)/BLKSIZE));
  123.     for(i=0;i<BLKSIZE_s;i++)window_s[i]=0.5*(1-cos(2.0*PI*(i-0.5)/BLKSIZE_s));
  124.   }
  125.   
  126.   /**********************************************************************
  127.    *  Delay signal by sync_flush samples                             *
  128.    **********************************************************************/
  129.   for ( j = 0; j < sync_flush; j++ ) /* for long window samples */
  130.     savebuf[j] = savebuf[j+flush];
  131.   
  132.   for ( j = sync_flush; j < syncsize; j++ )
  133.     savebuf[j] = buffer[j-sync_flush];
  134.   
  135.   
  136.   
  137.   /**********************************************************************
  138.    *  compute FFTs
  139.    **********************************************************************/
  140.   for ( j = 0; j < BLKSIZE; j++ ) 
  141.     wsamp_r[j] = window[j] * savebuf[j] * scalefac;
  142.   fft( wsamp_r, energy, ax, bx, 1024 );        /**long FFT**/
  143.   
  144. #ifdef HAVEGTK
  145.   /* delay by one granule.  Note we have to 'toggle' gr to do this */
  146.   if(gtkflag) {
  147.     for (j=0; j<HBLKSIZE ; j++) {
  148.       pinfo->energy[1-gr][chn][j]=energy_save[chn][j];
  149.       energy_save[chn][j]=energy[j];
  150.     }
  151.   }
  152. #endif
  153.   for ( sblock = 0; sblock < 3; sblock++ ) {
  154.     int shlen = 128;
  155.     int shoff = 2;
  156.     if (gpsycho) {
  157.       shlen = 192;
  158.       shoff = 1;
  159.     }
  160.     for ( j = 0, k = shlen * (shoff + sblock); j < 256; j++, k++ )
  161.       wsamp_rs[j] = window_s[j]* savebuf[k]*scalefac;
  162.     fft( wsamp_rs, &energy_s[sblock][0], &ax_s[sblock][0], &bx_s[sblock][0], 256 );
  163.   }
  164.  
  165.   /**********************************************************************/
  166.   /* compute side_energy / total_energy ratio */
  167.   /**********************************************************************/
  168.   *ms_ener_ratio=0;
  169.   if (check_ms_stereo) {
  170.     double ms_ener_side;
  171.     
  172.     /* compute energy in frequencies above MSFREQ */
  173.     for (ms_ener[chn]=0, j=MSFREQ; j<HBLKSIZE ; j++) {
  174.       ms_ener[chn] += energy[j];
  175.     }
  176.     if (chn == 1) {
  177.       /* compute energy in side channel */
  178.       ms_ener_side = fft_side( wsamp_r_int, MSFREQ);
  179.       
  180.       /* ener_side / total_energy */
  181.       if (ms_ener_side <= .01*(ms_ener[0]+ms_ener[1]))
  182.     ms_ener_side = .01;
  183.       else
  184.     *ms_ener_ratio = ms_ener_side/(ms_ener[0]+ms_ener[1]);
  185.     }
  186.   }
  187.  
  188.   if (force_ms) {
  189.     for (ms_ener[chn]=0, j=MSFREQ; j<HBLKSIZE ; j++)
  190.       ms_ener[chn] += energy[j];
  191.     if (chn == 1) 
  192.       if (ms_ener[1] <= .01*(ms_ener[0]+ms_ener[1]))
  193.     *ms_ener_ratio = .01;
  194.       else
  195.     *ms_ener_ratio = ms_ener[1]/(ms_ener[0]+ms_ener[1]);
  196.   }
  197.       
  198.  
  199. }
  200.  
  201.  
  202.  
  203.  
  204.  
  205.  
  206.  
  207.  
  208.  
  209. void L3psycho_anal( short int *buffer[2], int stereo,
  210.             int gr , layer * info,
  211.             double sfreq, 
  212.             int check_ms_stereo, double *ms_ener_ratio,
  213.             double ratio_d[2][21], double ratio_ds[2][12][3],
  214.             double percep_energy[2], int blocktype_d[2])
  215. {
  216.   static double pe[2]={0,0};
  217.   static double ms_ratio_s_old=0,ms_ratio_l_old=0;
  218.  
  219.   static double ratio[2][SBMAX_l];
  220.   static double ratio_s[2][SBMAX_s][3];
  221. #ifdef HAVEGTK
  222.   static double pe_save[2];
  223.   static double ers_save[2];
  224. #endif
  225.   static double thm_save[2][SBMAX_l];
  226.   static double en_save[2][SBMAX_l];
  227.   static double thm_s_save[2][SBMAX_s][3];
  228.   static double en_s_save[2][SBMAX_s][3];
  229.  
  230.  
  231.   int blocktype[2],uselongblock[2],chn;
  232.   unsigned int   b, i, j, k;
  233.   FLOAT          freq_mult, bval_lo;
  234.   double         temp1,temp2,temp3;
  235.   double ms_ratio_l=0,ms_ratio_s=0;
  236.   
  237.   /*         nint(); Layer III */
  238.   double   thr[CBANDS];
  239.   FLOAT ax[HBLKSIZE], bx[HBLKSIZE];
  240.   FLOAT energy[HBLKSIZE];
  241.   FLOAT energy_s[3][HBLKSIZE_s];
  242.   FLOAT ax_s[3][HBLKSIZE_s], bx_s[3][HBLKSIZE_s]; /* 256 samples not 129. */
  243.  
  244.   static float mld_l[SBMAX_l],mld_s[SBMAX_s];
  245.   
  246.  
  247.   
  248.  
  249. /* The static variables "r", "phi_sav", "new", "old" and "oldest" have    */
  250. /* to be remembered for the unpredictability measure.  For "r" and        */
  251. /* "phi_sav", the first index from the left is the channel select and     */
  252. /* the second index is the "age" of the data.                             */
  253.  static int     new = 0, old = 1, oldest = 0;
  254.  static int     init = 0,  sfreq_idx;
  255.  static double     cw[HBLKSIZE], eb[CBANDS];
  256.  static double     ctb[CBANDS];
  257.  static double    SNR_l[CBANDS], SNR_s[CBANDS_s];
  258.  static int    init_L3;
  259.  static double    minval[CBANDS],qthr_l[CBANDS],norm_l[CBANDS];
  260.  static double    qthr_s[CBANDS_s],norm_s[CBANDS_s];
  261.  static double    nb_1[2][CBANDS], nb_2[2][CBANDS];
  262.  
  263. /* Scale Factor Bands */
  264.  static int    cbw_l[SBMAX_l],bu_l[SBMAX_l],bo_l[SBMAX_l] ;
  265.  static int    cbw_s[SBMAX_s],bu_s[SBMAX_s],bo_s[SBMAX_s] ;
  266.  static double    w1_l[SBMAX_l], w2_l[SBMAX_l];
  267.  static double    w1_s[SBMAX_s], w2_s[SBMAX_s];
  268.  static double    en[SBMAX_l],   thm[SBMAX_l] ;
  269.  static int    blocktype_old[2] ;
  270.  int    sb,sblock;
  271.  static int    partition_l[HBLKSIZE],partition_s[HBLKSIZE_s];
  272.  
  273.  
  274. /* The following static variables are constants.                           */
  275.  
  276.  static FLOAT   crit_band[27] = {0,  100,  200, 300, 400, 510, 630,  770,
  277.                                920, 1080, 1270,1480,1720,2000,2320, 2700,
  278.                               3150, 3700, 4400,5300,6400,7700,9500,12000,
  279.                              15500,25000,30000};
  280.  
  281.  static FLOAT          *nb, *cb, *ecb, *bc;
  282.  
  283.  static FLOAT          *c, *fthr;
  284.  static F32            *snrtmp;
  285.  static    int    *numlines ;
  286.  static    int    *numlines2 ;
  287.  static int     *partition;
  288.  static FLOAT   *cbval, *rnorm;
  289.  static FLOAT   *absthr;
  290.  static double  *tmn;
  291.  static FCB     *s;
  292.  static FHBLK   *lthr;
  293.  static F2HBLK  *ax_sav, *bx_sav, *rx_sav;
  294.  
  295.  
  296.  assert( info->lay == 3 );
  297.  if(init==0){
  298. /* These dynamic memory allocations simulate "static" variables placed    */
  299. /* in the data space.  Each mem_alloc() call here occurs only once at     */
  300. /* initialization time.  The mem_free() function must not be called.      */
  301.  nb = (FLOAT *) mem_alloc(sizeof(FCB), "nb");
  302.  cb = (FLOAT *) mem_alloc(sizeof(FCB), "cb");
  303.  ecb = (FLOAT *) mem_alloc(sizeof(FCB), "ecb");
  304.  bc = (FLOAT *) mem_alloc(sizeof(FCB), "bc");
  305.  
  306.  c = (FLOAT *) mem_alloc(sizeof(FHBLK), "c");
  307.  fthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "fthr");
  308.  snrtmp = (F32 *) mem_alloc(sizeof(F2_32), "snrtmp");
  309.  
  310.      numlines = (int *) mem_alloc(sizeof(ICB), "numlines");
  311.      if (gpsycho) 
  312.        numlines2 = (int *) mem_alloc(sizeof(ICB), "numlines");
  313.      else  /* ISO model overwrites numlines with garbage */
  314.        numlines2 = numlines;
  315.  
  316.      partition = (int *) mem_alloc(sizeof(IHBLK), "partition");
  317.      cbval = (FLOAT *) mem_alloc(sizeof(FCB), "cbval");
  318.      rnorm = (FLOAT *) mem_alloc(sizeof(FCB), "rnorm");
  319.      absthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "absthr"); 
  320.      tmn = (double *) mem_alloc(sizeof(DCB), "tmn");
  321.      s = (FCB *) mem_alloc(sizeof(FCBCB), "s");
  322.      lthr = (FHBLK *) mem_alloc(sizeof(F2HBLK), "lthr");
  323.      rx_sav = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "rx_sav");
  324.      ax_sav = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "ax_sav");
  325.      bx_sav = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "bx_sav");
  326.  
  327.  
  328.      i = sfreq + 0.5;
  329.      switch(i){
  330.         case 32000: sfreq_idx = 0; break;
  331.         case 44100: sfreq_idx = 1; break;
  332.         case 48000: sfreq_idx = 2; break;
  333.         default:    printf("error, invalid sampling frequency: %d Hz\n",i);
  334.         exit(-1);
  335.      }
  336.      /* printf("absthr[][] sampling frequency index: %d\n",sfreq_idx); */
  337.      read_absthr(absthr, sfreq_idx);
  338.  
  339.  
  340. /* reset states used in unpredictability measure */
  341.     memset (rx_sav,0, 4* HBLKSIZE);
  342.     memset (ax_sav,0, 4* HBLKSIZE);
  343.     memset (bx_sav,0, 4* HBLKSIZE);
  344.      for(i=0;i<HBLKSIZE;i++){
  345.         lthr[0][i] = 60802371420160.0;
  346.         lthr[1][i] = 60802371420160.0;
  347.      }
  348. /*****************************************************************************
  349.  * Initialization: Compute the following constants for use later             *
  350.  *    partition[HBLKSIZE] = the partition number associated with each        *
  351.  *                          frequency line                                   *
  352.  *    cbval[CBANDS]       = the center (average) bark value of each          *
  353.  *                          partition                                        *
  354.  *    numlines[CBANDS]    = the number of frequency lines in each partition  *
  355.  *    tmn[CBANDS]         = tone masking noise                               *
  356.  *****************************************************************************/
  357. /* compute fft frequency multiplicand */
  358.      freq_mult = sfreq/BLKSIZE;
  359.  
  360. /* calculate fft frequency, then bval of each line (use fthr[] as tmp storage)*/
  361.      for(i=0;i<HBLKSIZE;i++){
  362.         temp1 = i*freq_mult;
  363.         j = 1;
  364.         while(temp1>crit_band[j])j++;
  365.         fthr[i]=j-1+(temp1-crit_band[j-1])/(crit_band[j]-crit_band[j-1]);
  366.      }
  367.      partition[0] = 0;
  368. /* temp2 is the counter of the number of frequency lines in each partition */
  369.      temp2 = 1;
  370.      cbval[0]=fthr[0];
  371.      bval_lo=fthr[0];
  372.      for(i=1;i<HBLKSIZE;i++){
  373.         if((fthr[i]-bval_lo)>0.33){
  374.            partition[i]=partition[i-1]+1;
  375.            cbval[partition[i-1]] = cbval[partition[i-1]]/temp2;
  376.            cbval[partition[i]] = fthr[i];
  377.            bval_lo = fthr[i];
  378.            numlines[partition[i-1]] = temp2;
  379.            temp2 = 1;
  380.         }
  381.         else {
  382.            partition[i]=partition[i-1];
  383.            cbval[partition[i]] += fthr[i];
  384.            temp2++;
  385.         }
  386.      }
  387.      numlines[partition[i-1]] = temp2;
  388.      cbval[partition[i-1]] = cbval[partition[i-1]]/temp2;
  389.  
  390. /************************************************************************
  391.  * Now compute the spreading function, s[j][i], the value of the spread-*
  392.  * ing function, centered at band j, for band i, store for later use    *
  393.  ************************************************************************/
  394.      for(j=0;j<CBANDS;j++){
  395.         for(i=0;i<CBANDS;i++){
  396.            temp1 = (cbval[i] - cbval[j])*1.05;
  397.            if(temp1>=0.5 && temp1<=2.5){
  398.               temp2 = temp1 - 0.5;
  399.               temp2 = 8.0 * (temp2*(temp2 - 2.0));
  400.            }
  401.            else temp2 = 0;
  402.            temp1 += 0.474;
  403.            temp3 = 15.811389+7.5*temp1-17.5*sqrt((double) (1.0+temp1*temp1));
  404.            if(temp3 <= -100) s[i][j] = 0;
  405.            else {
  406.               temp3 = (temp2 + temp3)*LN_TO_LOG10;
  407.               s[i][j] = exp(temp3);
  408.            }
  409.         }
  410.      }
  411.  
  412.   /* Calculate Tone Masking Noise values */
  413.      for(j=0;j<CBANDS;j++){
  414.         temp1 = 15.5 + cbval[j];
  415.         tmn[j] = (temp1>24.5) ? temp1 : 24.5;
  416.   /* Calculate normalization factors for the net spreading functions */
  417.         rnorm[j] = 0;
  418.         for(i=0;i<CBANDS;i++){
  419.            rnorm[j] += s[j][i];
  420.         }
  421.      }
  422.  
  423.  
  424.     /* setup stereo demasking thresholds */
  425.     /* formula reverse enginerred from plot in paper */
  426.     for ( sb = 0; sb < SBMAX_s; sb++ ) {
  427.       double mld = 1.25*(1-cos(3.14159*sb/SBMAX_s))-2.5;
  428.       mld_s[sb] = pow(10.0,mld);
  429.     }
  430.     for ( sb = 0; sb < SBMAX_l; sb++ ) {
  431.       double mld = 1.25*(1-cos(3.14159*sb/SBMAX_l))-2.5;
  432.       mld_l[sb] = pow(10.0,mld);
  433.     }
  434.   
  435.  
  436.  
  437.  
  438.      init++;
  439.  }
  440.  
  441.  if ( init_L3 == 0 )
  442.    {
  443.      L3para_read( sfreq,numlines2,partition_l,minval,qthr_l,norm_l,s3_l,
  444.           partition_s,qthr_s,norm_s,SNR_s,
  445.           cbw_l,bu_l,bo_l,w1_l,w2_l, cbw_s,bu_s,bo_s,w1_s,w2_s );
  446.      init_L3 ++ ;
  447.    }
  448. /************************* End of Initialization *****************************/
  449.  
  450.  
  451.  
  452.  
  453.  
  454. #ifdef AWS
  455.  /* reduce switch_pe if there where preecho events in previous granules */
  456.  if (gpsycho)
  457.  {
  458.    int prev_granule_used_shortblock=0;
  459.    
  460.    for (chn=0; chn<stereo; chn++)
  461.      if (blocktype_old[chn]==SHORT_TYPE) 
  462.        prev_granule_used_shortblock=1;
  463.    
  464.    if (prev_granule_used_shortblock)
  465.      switch_pe=maximum(switch_pe-700,900);
  466.    else  switch_pe = minimum(switch_pe+200,1800);
  467.  }
  468. #endif
  469.  
  470.  
  471.  
  472.  
  473.  
  474.  
  475.  
  476.  
  477.  
  478.  
  479.  
  480.   
  481.  for (chn=0; chn<stereo; chn++) {
  482.  
  483.  
  484.  for ( j = 0; j < 21; j++ )
  485.    ratio_d[chn][j] = ratio[chn][j];
  486.  for ( j = 0; j < 12; j++ )
  487.    for ( i = 0; i < 3; i++ )
  488.      ratio_ds[chn][j][i] = ratio_s[chn][j][i];
  489.  /* buggy ISO model doesn't delay pe */
  490.  if (gpsycho) percep_energy[chn] = pe[chn]; 
  491.  
  492.  
  493.  /*if ( chn == 0 )*/
  494.    if ( new == 0 )
  495.      {
  496.        new = 1;
  497.        old = 0;
  498.        oldest = 1;
  499.      }
  500.    else
  501.      {
  502.        new = 0;
  503.        old = 1;
  504.        oldest = 0;
  505.      }
  506.  
  507.  
  508. /**********************************************************************
  509. *  compute FFTs
  510. **********************************************************************/
  511.  L3psycho_energy( buffer[chn], energy, ax, bx, energy_s, ax_s, bx_s,
  512.           chn,gr,info,check_ms_stereo,&ms_ratio_l);
  513.  ms_ratio_s=ms_ratio_l;
  514.  
  515. /**********************************************************************
  516. *    compute unpredicatability of first six spectral lines            * 
  517. **********************************************************************/
  518.  for ( j = 0; j < 6; j++ )
  519.    {     /* calculate unpredictability measure cw */
  520.      double an, a1, a2;
  521.      double bn, b1, b2;
  522.      double rn, r1, r2;
  523.      double numre, numim, den;
  524.  
  525.      a2 = ax_sav[chn][oldest][j];
  526.      b2 = bx_sav[chn][oldest][j];
  527.      r2 = rx_sav[chn][oldest][j];
  528.      a1 = ax_sav[chn][oldest][j] = ax_sav[chn][old][j];
  529.      b1 = bx_sav[chn][oldest][j] = bx_sav[chn][old][j];
  530.      r1 = rx_sav[chn][oldest][j] = rx_sav[chn][old][j];
  531.      an = ax_sav[chn][old][j] = ax[j];
  532.      bn = bx_sav[chn][old][j] = bx[j];
  533.      rn = rx_sav[chn][old][j] = sqrt(energy[j]);
  534.  
  535.      { /* square (x1,y1) */
  536.        if( r1 != 0.0 ) {
  537.      numre = (a1*b1);
  538.      numim = (a1*a1-b1*b1)*0.5;
  539.      den = r1*r1;
  540.        } else {
  541.      numre = 1.0;
  542.      numim = 0.0;
  543.      den = 1.0;
  544.        }
  545.      }
  546.  
  547.      { /* multiply by (x2,-y2) */
  548.        if( r2 != 0.0 ) {
  549.      double tmp2 = (numim+numre)*(a2+b2)*0.5;
  550.      double tmp1 = -a2*numre+tmp2;
  551.      numre =       -b2*numim+tmp2;
  552.      numim = tmp1;
  553.      den *= r2;
  554.        } else {
  555.      /* do nothing */
  556.        }
  557.      }
  558.  
  559.      { /* r-prime factor */
  560.        double tmp = (2.0*r1-r2)/den;
  561.        numre *= tmp;
  562.        numim *= tmp;
  563.      }
  564.  
  565.      if( (den=rn+fabs(2.0*r1-r2)) != 0.0 ) {
  566.        numre = (an+bn)/2.0-numre;
  567.        numim = (an-bn)/2.0-numim;
  568.        cw[j] = sqrt(numre*numre+numim*numim)/den;
  569.      } else {
  570.        cw[j] = 0.0;
  571.      }
  572.  
  573.    }
  574.                                                                                   
  575. /**********************************************************************
  576. *     compute unpredicatibility of next 200 spectral lines            *
  577. **********************************************************************/ 
  578.  sblock = 1;
  579.     
  580.  for ( j = 6; j < 206; j += 4 )
  581.    {/* calculate unpredictability measure cw */
  582.      double rn, r1, r2;
  583.      double numre, numim, den;
  584.  
  585.      k = (j+2) / 4; 
  586.  
  587.      { /* square (x1,y1) */
  588.        r1 = sqrt((double)energy_s[0][k]);
  589.        if( r1 != 0.0 ) {
  590.      double a1 = ax_s[0][k];
  591.      double b1 = bx_s[0][k];
  592.      numre = (a1*b1);
  593.      numim = (a1*a1-b1*b1)*0.5;
  594.      den = r1*r1;
  595.        } else {
  596.      numre = 1.0;
  597.      numim = 0.0;
  598.      den = 1.0;
  599.        }
  600.      }
  601.  
  602.      { /* multiply by (x2,-y2) */
  603.        r2 = sqrt((double)energy_s[2][k]);
  604.        if( r2 != 0.0 ) {
  605.      double a2 = ax_s[2][k];
  606.      double b2 = bx_s[2][k];
  607.  
  608.      double tmp2 = (numim+numre)*(a2+b2)*0.5;
  609.      double tmp1 = -a2*numre+tmp2;
  610.      numre =       -b2*numim+tmp2;
  611.      numim = tmp1;
  612.  
  613.      den *= r2;
  614.        } else {
  615.      /* do nothing */
  616.        }
  617.      }
  618.  
  619.      { /* r-prime factor */
  620.        double tmp = (2.0*r1-r2)/den;
  621.        numre *= tmp;
  622.        numim *= tmp;
  623.      }
  624.  
  625.      rn = sqrt((double)energy_s[1][k]);
  626.      if( (den=rn+fabs(2.0*r1-r2)) != 0.0 ) {
  627.        double an = ax_s[1][k];
  628.        double bn = bx_s[1][k];
  629.        numre = (an+bn)/2.0-numre;
  630.        numim = (an-bn)/2.0-numim;
  631.        cw[j] = sqrt(numre*numre+numim*numim)/den;
  632.      } else {
  633.        cw[j] = 0.0;
  634.      }
  635.  
  636.      cw[j+1] = cw[j+2] = cw[j+3] = cw[j];
  637.    }
  638.  
  639.  
  640. /**********************************************************************
  641. *    Set unpredicatiblility of remaining spectral lines to 0.4  206..513 *
  642. **********************************************************************/
  643.  for ( j = 206; j < HBLKSIZE; j++ )
  644.    cw[j] = 0.4;
  645.     
  646.  
  647.  
  648. #if 0
  649.  for ( j = 14; j < HBLKSIZE-4; j += 4 )
  650.    {/* calculate energy from short ffts */
  651.      double tot,ave;
  652.      k = (j+2) / 4; 
  653.      for (tot=0, sblock=0; sblock < 3; sblock++)
  654.        tot+=energy_s[sblock][k];
  655.      ave = energy[j+1]+ energy[j+2]+ energy[j+3]+ energy[j];
  656.      ave /= 4.;
  657.      /*
  658.        printf("energy / tot %i %5.2f   %e  %e\n",j,ave/(tot*16./3.),
  659.        ave,tot*16./3.);
  660.      */
  661.      energy[j+1] = energy[j+2] = energy[j+3] =  energy[j]=tot;
  662.    }
  663. #endif
  664.  
  665.  
  666.  
  667.  
  668.  
  669.  
  670.  
  671.  
  672. /**********************************************************************
  673. *    Calculate the energy and the unpredictability in the threshold   *
  674. *    calculation partitions                                           *
  675. **********************************************************************/
  676.  for ( b = 0; b < CBANDS; b++ )
  677.    {
  678.      eb[b] = 0.0;
  679.      cb[b] = 0.0;
  680.    }
  681.  for ( j = 0; j < HBLKSIZE; j++ )
  682.    {
  683.      int tp = partition_l[j];
  684.      if ( tp >= 0 )
  685.        {
  686.      eb[tp] += energy[j];
  687.      cb[tp] += cw[j] * energy[j];
  688.        }
  689.    }
  690.  
  691.  
  692. /**********************************************************************
  693. *      convolve the partitioned energy and unpredictability           *
  694. *      with the spreading function, s3_l[b][k]                        *
  695. ******************************************************************** */
  696.  for ( b = 0; b < CBANDS; b++ )
  697.    {
  698.      ecb[b] = 0.0;
  699.      ctb[b] = 0.0;
  700.    }
  701.  /* s3_l is a sparse matrix for 44.1khz */
  702.  if (sfreq_idx==1) {
  703.    sprdngf1(ecb,eb);
  704.    sprdngf2(ctb,cb);
  705.  }
  706.  else {
  707.    for ( b = 0;b < CBANDS; b++ )
  708.      {
  709.        for ( k = 0; k < CBANDS; k++ )
  710.      {
  711.        ecb[b] += s3_l[b][k] * eb[k];    /* sprdngf for Layer III */
  712.        ctb[b] += s3_l[b][k] * cb[k];
  713.      }
  714.      }
  715.  }
  716.  /* calculate the tonality of each threshold calculation partition */
  717.  /* calculate the SNR in each threshhold calculation partition */
  718.  
  719.  for ( b = 0; b < CBANDS; b++ )
  720.    {
  721.      double cbb,tbb;
  722.      if (ecb[b] != 0.0 )
  723.        {
  724.      cbb = ctb[b]/ecb[b];
  725.      if (cbb <0.01) cbb = 0.01;
  726.      cbb = log( cbb);
  727.        }
  728.      else
  729.        cbb = 0.0 ;
  730.      tbb = -0.299 - 0.43*cbb;  /* conv1=-0.299, conv2=-0.43 */
  731.      tbb = minimum( 1.0, maximum( 0.0, tbb) ) ;  /* 0<tbb<1 */
  732.      SNR_l[b] = maximum( minval[b], 29.0*tbb+6.0*(1.0-tbb) );
  733.    }    /* TMN=29.0,NMT=6.0 for all calculation partitions */
  734.  
  735.  for ( b = 0; b < CBANDS; b++ ) /* calculate the threshold for each partition */
  736.    nb[b] = ecb[b] * norm_l[b] * exp( -SNR_l[b] * LN_TO_LOG10 );
  737.  
  738.  for ( b = 0; b < CBANDS; b++ )
  739.    { /* pre-echo control */
  740.      double temp_1; /* BUG of IS */
  741.      int rpelev=2; int rpelev2=16; 
  742.      temp_1 = minimum( nb[b], minimum(rpelev*nb_1[chn][b],rpelev2*nb_2[chn][b]) );
  743.      thr[b] = maximum( qthr_l[b], temp_1 );/* rpelev=2.0, rpelev2=16.0 */
  744.  
  745.      nb_2[chn][b] = nb_1[chn][b];
  746.      nb_1[chn][b] = nb[b];
  747.    }
  748.  /* note: all surges in PE are because of the above pre-echo formula
  749.   * for temp_1.  it this is not used, PE is always around 600
  750.   */
  751.  
  752.  pe[chn] = 0.0;        /*  calculate percetual entropy */
  753.  for ( b = 0; b < CBANDS; b++ )
  754.    {
  755.      double tp ;
  756.      tp = minimum( 0.0, log((thr[b]+1.0) / (eb[b]+1.0) ) ) ;  /*not log*/
  757.      pe[chn] -= numlines[b] * tp ;
  758.  
  759.    }    /* thr[b] -> thr[b]+1.0 : for non sound portition */
  760.  if (!gpsycho) percep_energy[chn]=pe[chn];
  761.  
  762.  
  763.  /*************************************************************** 
  764.   * Check to see if we also need to compute long block thresholds
  765.   ***************************************************************/
  766.  uselongblock[chn] = (pe[chn] < switch_pe);
  767.  
  768.  {
  769.    double tot[3]={0,0,0};
  770.    double mn,mx;
  771.    for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
  772.      for (sblock=0; sblock < 3; sblock++)
  773.        tot[sblock]+=energy_s[sblock][j];
  774.    mn = minimum(tot[0],tot[1]);
  775.    mn = minimum(mn,tot[2]);
  776.    mx = maximum(tot[0],tot[1]);
  777.    mx = maximum(mx,tot[2]);
  778. #ifdef HAVEGTK
  779.    /* delay by one granule.  Note we have to 'toggle' gr to do this */
  780.    if (gtkflag) {
  781.      pinfo->ers[1-gr][chn]=ers_save[chn];
  782.      ers_save[chn]=mx/(1e-12+mn);
  783.      pinfo->pe[1-gr][chn]=pe_save[chn];
  784.      pe_save[chn]=pe[chn];
  785.    }
  786. #endif
  787. #ifdef ENER_AWS
  788.    if (gpsycho) {
  789.    uselongblock[chn] = 1;
  790.  
  791.    /* big surge of energy - always use short blocks */
  792.    if (  mx > 25*mn) uselongblock[chn] = 0;
  793.  
  794.    /* switch based on pe, but only if there is some energy surge */
  795.    if (( mx > 2.5*mn ) && (pe[chn] > switch_pe)) uselongblock[chn]=0;
  796.    }
  797. #endif
  798.  }
  799.  
  800.  
  801.  
  802.  
  803.  
  804.  if ((chn == 1) && (force_ms)) {
  805.    /* Forced ms_stereo mode.  */
  806.    /* ch=0 (mid) blocktype determines ch=1 (side) blocktype */
  807.    uselongblock[1] = uselongblock[0];
  808.  }
  809.  
  810.  
  811.  /*************************************************************** 
  812.   * compute masking thresholds
  813.   * we need to always compute short block masking thresholds because
  814.   * we dont know if this granule will be changed to a short block later.
  815.   ***************************************************************/
  816.  if ( uselongblock[chn]) 
  817.    {                /* no attack : use long blocks */
  818.      /* threshold calculation (part 2) */
  819.      for ( sb = 0; sb < SBMAX_l; sb++ )
  820.        {
  821.      en[sb] = w1_l[sb] * eb[bu_l[sb]] + w2_l[sb] * eb[bo_l[sb]];
  822.      thm[sb] = w1_l[sb] *thr[bu_l[sb]] + w2_l[sb] * thr[bo_l[sb]];
  823.      for ( b = bu_l[sb]+1; b < bo_l[sb]; b++ )
  824.        {
  825.          en[sb]  += eb[b];
  826.          thm[sb] += thr[b];
  827.        }
  828.      if ( en[sb] != 0.0 )
  829.        ratio[chn][sb] = thm[sb]/en[sb];
  830.      else
  831.        ratio[chn][sb] = 0.0;
  832.        }
  833. #ifdef HAVEGTK
  834.   /* delay by one granule.  Note we have to 'toggle' gr to do this */
  835.      if (gtkflag) {
  836.        for (sb=0; sb< SBMAX_l; sb ++ ) {
  837.      pinfo->thr[1-gr][chn][sb]=thm_save[chn][sb];
  838.      pinfo->en[1-gr][chn][sb]=en_save[chn][sb];
  839.      thm_save[chn][sb]=thm[sb];
  840.      en_save[chn][sb]=en[sb];
  841.        }
  842.      }
  843. #endif
  844.      for (sb=0; sb< SBMAX_l; sb ++ ) {
  845.        thm_save[chn][sb]=thm[sb];
  846.        en_save[chn][sb]=en[sb];
  847.      }
  848.  
  849.    }
  850.  
  851.  /* in some cases, when computing the next granule, we may switch this
  852.   * granule to a short block.  compute short block thresholds just in case */
  853.  /*else*/
  854.  if ( gpsycho  || !uselongblock[chn])
  855.    {
  856.      
  857.      /* threshold calculation for short blocks */
  858.      
  859.      for ( sblock = 0; sblock < 3; sblock++ )
  860.        {
  861.      for ( b = 0; b < CBANDS_s; b++ )
  862.        {
  863.          eb[b] = 0.0;
  864.          ecb[b] = 0.0;
  865.        }
  866.      for ( j = 0; j < HBLKSIZE_s; j++ )
  867.        eb[partition_s[j]] += energy_s[sblock][j];
  868.      for ( b = 0; b < CBANDS_s; b++ )
  869.        for ( k = 0; k < CBANDS_s; k++ )
  870.          ecb[b] += s3_l[b][k] * eb[k];
  871.      for ( b = 0; b < CBANDS_s; b++ )
  872.        {
  873.          if (gpsycho)
  874.            nb[b] = ecb[b] * norm_s[b] * exp( (double) SNR_s[b] * LN_TO_LOG10 );
  875.          else
  876.            nb[b] = ecb[b] * norm_l[b] * exp( (double) SNR_s[b] * LN_TO_LOG10 );
  877.          thr[b] = maximum (qthr_s[b],nb[b]);
  878.        }
  879.      for ( sb = 0; sb < SBMAX_s; sb++ )
  880.        {
  881.          en[sb] = w1_s[sb] * eb[bu_s[sb]] + w2_s[sb] * eb[bo_s[sb]];
  882.          thm[sb] = w1_s[sb] *thr[bu_s[sb]] + w2_s[sb] * thr[bo_s[sb]];
  883.          for ( b = bu_s[sb]+1; b < bo_s[sb]; b++ )
  884.            {
  885.          en[sb] += eb[b];
  886.          thm[sb] += thr[b];
  887.            }
  888.          if ( en[sb] != 0.0 )
  889.            ratio_s[chn][sb][sblock] = thm[sb]/en[sb];
  890.          else
  891.            ratio_s[chn][sb][sblock] = 0.0;
  892. #ifdef HAVEGTK
  893.   /* delay by one granule.  Note we have to 'toggle' gr to do this */
  894.          if (gtkflag) {
  895.            pinfo->thr_s[1-gr][chn][3*sb+sblock]=thm_s_save[chn][sb][sblock];
  896.            pinfo->en_s[1-gr][chn][3*sb+sblock]=en_s_save[chn][sb][sblock];
  897.          }
  898. #endif
  899.          thm_s_save[chn][sb][sblock]=thm[sb];
  900.          en_s_save[chn][sb][sblock]=en[sb];
  901.        }
  902.        }
  903.    } 
  904.  
  905.  
  906.  /* compute M/S thresholds from Johnston & Ferreira 1992 ICASSP paper */
  907.  if (force_ms) {
  908.  
  909.    if (chn == 1) {
  910.      double rside,rmid,mld; 
  911.      if ( uselongblock[chn]) {
  912.        for ( sb = 0; sb < SBMAX_l; sb++ ) {
  913.      mld = mld_l[sb];
  914.      rmid = Max(ratio[0][sb],Min(ratio[1][sb],mld));
  915.      rside = Max(ratio[1][sb],Min(ratio[0][sb],mld));
  916.      ratio[0][sb]=rmid;
  917.      ratio[1][sb]=rside;
  918.        }
  919.      } 
  920.      /* alwasy compute these - we may need them later */
  921.      for ( sblock = 0; sblock < 3; sblock++ ){
  922.        for ( sb = 0; sb < SBMAX_s; sb++ ) {
  923.      mld = mld_s[sb];
  924.      rmid = Max(ratio_s[0][sb][sblock],Min(ratio_s[1][sb][sblock],mld));
  925.      rside = Max(ratio_s[1][sb][sblock],Min(ratio_s[0][sb][sblock],mld));
  926.      ratio_s[0][sb][sblock]=rmid;
  927.      ratio_s[1][sb][sblock]=rside;
  928.        }
  929.      }
  930.    }
  931.  }
  932.  } /* end loop over chn */
  933.  
  934.  
  935.  
  936. #if 0
  937.  /* determin ms_ratio from masking thresholds*/
  938.  { double sidetot=0,tot=0,side,mid;
  939.  for (sb= 15 ; sb< SBMAX_l; sb ++ ) {
  940.    side = thm_save[0][sb]-thm_save[1][sb];
  941.    mid = thm_save[0][sb]+thm_save[1][sb];
  942.    sidetot += side*side;
  943.    tot += side*side + mid*mid;
  944.  }
  945.  ms_ratio_l = Max(.875*ms_ratio_l,1.2*sidetot/(1e-10+tot));
  946.  
  947.  sidetot=0; tot=0;
  948.  for ( sblock = 0; sblock < 3; sblock++ )
  949.    for ( sb = 9; sb < SBMAX_s; sb++ ) {
  950.      side = thm_s_save[0][sb][sblock]-thm_s_save[1][sb][sblock];
  951.      mid = thm_s_save[0][sb][sblock]+thm_s_save[1][sb][sblock];
  952.      sidetot += side*side;
  953.      tot += side*side + mid*mid;
  954.    }
  955.  ms_ratio_s = Max(.875*ms_ratio_s,1.2*sidetot/(1e-10+tot));
  956.  }
  957. #endif
  958.  
  959.  
  960.  
  961.  /*************************************************************** 
  962.   * determin final block type
  963.   ***************************************************************/
  964.  for (chn=0; chn<stereo; chn++) {
  965.    blocktype[chn] = NORM_TYPE;
  966.  }
  967.  
  968.  if (!allow_diff_short)
  969.  if (info->mode==MPG_MD_JOINT_STEREO) {
  970.    /* force both channels to use the same block type */
  971.    /* this is necessary if the frame is to be encoded in ms_stereo.  */
  972.    /* But even without ms_stereo, FhG  does this */
  973.    int bothlong= (uselongblock[0] && uselongblock[1]);
  974.    if (!bothlong) {
  975.      uselongblock[0]=0;
  976.      uselongblock[1]=0;
  977.    }
  978.  }
  979.  
  980.  
  981.  
  982.  
  983.  /* update the blocktype of the previous granule, since it depends on what
  984.   * happend in this granule */
  985.  for (chn=0; chn<stereo; chn++) {
  986.  if ( uselongblock[chn])
  987.    {                /* no attack : use long blocks */
  988.      switch( blocktype_old[chn] ) 
  989.        {
  990.        case NORM_TYPE:
  991.        case STOP_TYPE:
  992.      blocktype[chn] = NORM_TYPE;
  993.      break;
  994.        case SHORT_TYPE:
  995.      blocktype[chn] = STOP_TYPE; 
  996.      break;
  997.        case START_TYPE:
  998.      fprintf( stderr, "Error in block selecting\n" );
  999.      abort();
  1000.      break; /* problem */
  1001.        }
  1002.    } else   {
  1003.      /* attack : use short blocks */
  1004.      blocktype[chn] = SHORT_TYPE;
  1005.      if ( blocktype_old[chn] == NORM_TYPE ) {
  1006.        blocktype_old[chn] = START_TYPE;
  1007.      }
  1008.      if ( blocktype_old[chn] == STOP_TYPE ) {
  1009.        blocktype_old[chn] = SHORT_TYPE ;
  1010.      }
  1011.    }
  1012.  
  1013.  
  1014.  blocktype_d[chn] = blocktype_old[chn];  /* value returned to calling program */
  1015.  blocktype_old[chn] = blocktype[chn];    /* save for next call to l3psy_anal */
  1016.  }
  1017.  
  1018.  if (blocktype_d[0]==2)
  1019.    *ms_ener_ratio = ms_ratio_s_old;
  1020.  else
  1021.    *ms_ener_ratio = ms_ratio_l_old;
  1022.  
  1023.  ms_ratio_s_old = ms_ratio_s;
  1024.  ms_ratio_l_old = ms_ratio_l;
  1025. }
  1026.  
  1027.  
  1028.  
  1029.  
  1030.  
  1031.  
  1032. extern double psy_data[];
  1033.  
  1034. void L3para_read(double sfreq, int *numlines, int *partition_l, double *minval, double *qthr_l, double *norm_l, double (*s3_l)[63], int *partition_s, double *qthr_s, double *norm_s, double *SNR, int *cbw_l, int *bu_l, int *bo_l, double *w1_l, double *w2_l, int *cbw_s, int *bu_s, int *bo_s, double *w1_s, double *w2_s)
  1035. {
  1036.   double freq_tp;
  1037.   static double bval_l[CBANDS], bval_s[CBANDS];
  1038.   int   cbmax=0, cbmax_tp;
  1039.   static double s3_s[CBANDS][CBANDS];
  1040.   double *p = psy_data;
  1041.  
  1042.   int  sbmax ;
  1043.   int  i,j,k,k2,loop, part_max ;
  1044.  
  1045.   /* Read long block data */
  1046.  
  1047.   for(loop=0;loop<6;loop++)
  1048.     {
  1049.       freq_tp = *p++;
  1050.       cbmax_tp = (int) *p++;
  1051.       cbmax_tp++;
  1052.  
  1053.       if (sfreq == freq_tp )
  1054.     {
  1055.       cbmax = cbmax_tp;
  1056.       for(i=0,k2=0;i<cbmax_tp;i++)
  1057.         {
  1058.           j = (int) *p++;
  1059.           numlines[i] = (int) *p++;
  1060.           minval[i] = *p++;
  1061.           qthr_l[i] = *p++;
  1062.           norm_l[i] = *p++;
  1063.           bval_l[i] = *p++;
  1064.           if (j!=i)
  1065.         {
  1066.           printf("1. please check \"psy_data\"");
  1067.           exit(-1);
  1068.         }
  1069.           for(k=0;k<numlines[i];k++)
  1070.         partition_l[k2++] = i ;
  1071.         }
  1072.     }
  1073.       else
  1074.     p += cbmax_tp * 6;
  1075.     }
  1076.  
  1077.   /************************************************************************
  1078.    * Now compute the spreading function, s[j][i], the value of the spread-*
  1079.    * ing function, centered at band j, for band i, store for later use    *
  1080.    ************************************************************************/
  1081.   part_max = cbmax ;
  1082.   for(i=0;i<part_max;i++)
  1083.     {
  1084.       double tempx,x,tempy,temp;
  1085.       for(j=0;j<part_max;j++)
  1086.     {
  1087.       tempx = (bval_l[i] - bval_l[j])*1.05;
  1088.       if (j>=i) tempx = (bval_l[i] - bval_l[j])*3.0;
  1089.       else    tempx = (bval_l[i] - bval_l[j])*1.5;
  1090.       /*             if (j>=i) tempx = (bval_l[j] - bval_l[i])*3.0;
  1091.              else    tempx = (bval_l[j] - bval_l[i])*1.5; */
  1092.       if(tempx>=0.5 && tempx<=2.5)
  1093.         {
  1094.           temp = tempx - 0.5;
  1095.           x = 8.0 * (temp*temp - 2.0 * temp);
  1096.         }
  1097.       else x = 0.0;
  1098.       tempx += 0.474;
  1099.       tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
  1100.       if (tempy <= -60.0) s3_l[i][j] = 0.0;
  1101.       else                s3_l[i][j] = exp( (x + tempy)*LN_TO_LOG10 );
  1102.     }
  1103.     }
  1104.  
  1105.   /* Read short block data */
  1106.  
  1107.   for(loop=0;loop<6;loop++)
  1108.     {
  1109.       freq_tp = *p++;
  1110.       cbmax_tp = (int) *p++;
  1111.       cbmax_tp++;
  1112.  
  1113.       if (sfreq == freq_tp )
  1114.     {
  1115.       cbmax = cbmax_tp;
  1116.       for(i=0,k2=0;i<cbmax_tp;i++)
  1117.         {
  1118.           j = (int) *p++;
  1119.           numlines[i] = (int) *p++;
  1120.           qthr_s[i] = *p++;
  1121.           norm_s[i] = *p++;
  1122.           SNR[i] = *p++;
  1123.           bval_s[i] = *p++;
  1124.           if (j!=i)
  1125.         {
  1126.           printf("3. please check \"psy_data\"");
  1127.           exit(-1);
  1128.         }
  1129.           for(k=0;k<numlines[i];k++)
  1130.         partition_s[k2++] = i ;
  1131.         }
  1132.     }
  1133.       else
  1134.     p += cbmax_tp * 6;
  1135.     }
  1136.  
  1137.   /************************************************************************
  1138.    * Now compute the spreading function, s[j][i], the value of the spread-*
  1139.    * ing function, centered at band j, for band i, store for later use    *
  1140.    ************************************************************************/
  1141.   part_max = cbmax ;
  1142.   for(i=0;i<part_max;i++)
  1143.     {
  1144.       double tempx,x,tempy,temp;
  1145.       for(j=0;j<part_max;j++)
  1146.     {
  1147.       tempx = (bval_s[i] - bval_s[j])*1.05;
  1148.       if (j>=i) tempx = (bval_s[i] - bval_s[j])*3.0;
  1149.       else    tempx = (bval_s[i] - bval_s[j])*1.5;
  1150.       if(tempx>=0.5 && tempx<=2.5)
  1151.         {
  1152.           temp = tempx - 0.5;
  1153.           x = 8.0 * (temp*temp - 2.0 * temp);
  1154.         }
  1155.       else x = 0.0;
  1156.       tempx += 0.474;
  1157.       tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
  1158.       if (tempy <= -60.0) s3_s[i][j] = 0.0;
  1159.       else                s3_s[i][j] = exp( (x + tempy)*LN_TO_LOG10 );
  1160.     }
  1161.     }
  1162.   /* Read long block data for converting threshold calculation 
  1163.      partitions to scale factor bands */
  1164.  
  1165.   for(loop=0;loop<6;loop++)
  1166.     {
  1167.       freq_tp = *p++;
  1168.       sbmax =  (int) *p++;
  1169.       sbmax++;
  1170.  
  1171.       if (sfreq == freq_tp)
  1172.     {
  1173.       for(i=0;i<sbmax;i++)
  1174.         {
  1175.           j = (int) *p++;
  1176.           cbw_l[i] = (int) *p++;
  1177.           bu_l[i] = (int) *p++;
  1178.           bo_l[i] = (int) *p++;
  1179.           w1_l[i] = (double) *p++;
  1180.           w2_l[i] = (double) *p++;
  1181.           if (j!=i)
  1182.         { printf("30:please check \"psy_data\"\n");
  1183.         exit(-1);
  1184.         }
  1185.  
  1186.           if (i!=0)
  1187.         if ( (bo_l[i] != (bu_l[i]+cbw_l[i])) ||
  1188.              (fabs(1.0-w1_l[i]-w2_l[i-1]) > 0.01 ) )
  1189.           {
  1190.             printf("31l: please check \"psy_data.\"\n");
  1191.             exit(-1);
  1192.           }
  1193.         }
  1194.     }
  1195.       else
  1196.     p += sbmax * 6;
  1197.     }
  1198.  
  1199.   /* Read short block data for converting threshold calculation 
  1200.      partitions to scale factor bands */
  1201.  
  1202.   for(loop=0;loop<6;loop++)
  1203.     {
  1204.       freq_tp = *p++;
  1205.       sbmax = (int) *p++;
  1206.       sbmax++;
  1207.  
  1208.       if (sfreq == freq_tp)
  1209.     {
  1210.       for(i=0;i<sbmax;i++)
  1211.         {
  1212.           j = (int) *p++;
  1213.           cbw_s[i] = (int) *p++;
  1214.           bu_s[i] = (int) *p++;
  1215.           bo_s[i] = (int) *p++;
  1216.           w1_s[i] = *p++;
  1217.           w2_s[i] = *p++;
  1218.           if (j!=i)
  1219.         { printf("30:please check \"psy_data\"\n");
  1220.         exit(-1);
  1221.         }
  1222.  
  1223.           if (i!=0)
  1224.         if ( (bo_s[i] != (bu_s[i]+cbw_s[i])) ||
  1225.              (fabs(1.0-w1_s[i]-w2_s[i-1]) > 0.01 ) )
  1226.           { printf("31s: please check \"psy_data.\"\n");
  1227.           exit(-1);
  1228.           }
  1229.         }
  1230.     }
  1231.       else
  1232.     p += sbmax * 6;
  1233.     }
  1234.  
  1235. }
  1236. void sprdngf1(dest,source)
  1237. FLOAT *dest;
  1238. double *source;
  1239. {
  1240. #ifdef GOOBER
  1241. static int mfcinit = 0;
  1242. #endif
  1243. int b,k;
  1244. extern double s3_l[CBANDS][CBANDS];  /* s3_l is a sparse matrix. fudge with it for 44.1 khz */
  1245.  static int s3ind[CBANDS][2] = {
  1246.    {0,2},
  1247.    {0,3},
  1248.    {0,4},
  1249.    {0,5},
  1250.    {0,6},
  1251.    {0,7},
  1252.    {0,8},
  1253.    {0,9},
  1254.    {0,10},
  1255.    {0,11},
  1256.    {0,12},
  1257.    {1,14},
  1258.    {1,14},
  1259.    {2,15},
  1260.    {3,15},
  1261.    {5,16},
  1262.    {6,17},
  1263.    {7,19},
  1264.    {9,20},
  1265.    {10,21},
  1266.    {11,22},
  1267.    {12,23},
  1268.    {14,24},
  1269.    {15,25},
  1270.    {15,27},
  1271.    {16,28},
  1272.    {16,28},
  1273.    {17,29},
  1274.    {18,30},
  1275.    {19,31},
  1276.    {19,32},
  1277.    {20,34},
  1278.    {21,35},
  1279.    {22,36},
  1280.    {22,36},
  1281.    {23,37},
  1282.    {24,38},
  1283.    {25,39},
  1284.    {26,41},
  1285.    {27,42},
  1286.    {28,43},
  1287.    {29,44},
  1288.    {30,45},
  1289.    {31,46},
  1290.    {32,47},
  1291.    {33,48},
  1292.    {34,49},
  1293.    {35,50},
  1294.    {36,51},
  1295.    {37,52},
  1296.    {37,53},
  1297.    {38,54},
  1298.    {39,55},
  1299.    {40,56},
  1300.    {41,57},
  1301.    {42,58},
  1302.    {43,59},
  1303.    {44,60},
  1304.    {45,61},
  1305.    {46,62},
  1306.    {47,62},
  1307.    {48,62},
  1308.    {48,62},
  1309.  };
  1310.  
  1311.  
  1312.     for ( b = 0;b < CBANDS; b++ )
  1313.         for ( k = s3ind[b][0]; k <= s3ind[b][1]; k++ )
  1314.         dest[b] += s3_l[b][k] * source[k];    
  1315. #ifdef GOOBER
  1316.     if (mfcinit==0) {
  1317.       mfcinit++;
  1318.     for (b=0;b<CBANDS;b++) {
  1319.       fprintf("{",b);
  1320.       for (k=0;k<CBANDS;k++) {
  1321.         if (fabs(s3_l[b][k])>0.00000001)
  1322.           fprintf("%i ",k);
  1323.       }
  1324.       fprintf("\n");
  1325.     }
  1326.     } 
  1327. #endif /* GOOBER */
  1328.  
  1329. }
  1330.  
  1331.  
  1332. void sprdngf2(dest,source)
  1333. double *dest;
  1334. FLOAT *source;
  1335. {
  1336. int b,k;
  1337. extern double s3_l[CBANDS][CBANDS];  /* s3_l is a sparse matrix. fudge with it for 44.1 khz */
  1338.  static int s3ind[CBANDS][2] = {
  1339.    {0,2},
  1340.    {0,3},
  1341.    {0,4},
  1342.    {0,5},
  1343.    {0,6},
  1344.    {0,7},
  1345.    {0,8},
  1346.    {0,9},
  1347.    {0,10},
  1348.    {0,11},
  1349.    {0,12},
  1350.    {1,14},
  1351.    {1,14},
  1352.    {2,15},
  1353.    {3,15},
  1354.    {5,16},
  1355.    {6,17},
  1356.    {7,19},
  1357.    {9,20},
  1358.    {10,21},
  1359.    {11,22},
  1360.    {12,23},
  1361.    {14,24},
  1362.    {15,25},
  1363.    {15,27},
  1364.    {16,28},
  1365.    {16,28},
  1366.    {17,29},
  1367.    {18,30},
  1368.    {19,31},
  1369.    {19,32},
  1370.    {20,34},
  1371.    {21,35},
  1372.    {22,36},
  1373.    {22,36},
  1374.    {23,37},
  1375.    {24,38},
  1376.    {25,39},
  1377.    {26,41},
  1378.    {27,42},
  1379.    {28,43},
  1380.    {29,44},
  1381.    {30,45},
  1382.    {31,46},
  1383.    {32,47},
  1384.    {33,48},
  1385.    {34,49},
  1386.    {35,50},
  1387.    {36,51},
  1388.    {37,52},
  1389.    {37,53},
  1390.    {38,54},
  1391.    {39,55},
  1392.    {40,56},
  1393.    {41,57},
  1394.    {42,58},
  1395.    {43,59},
  1396.    {44,60},
  1397.    {45,61},
  1398.    {46,62},
  1399.    {47,62},
  1400.    {48,62},
  1401.    {48,62},
  1402.  };
  1403.  
  1404.  
  1405.     for ( b = 0;b < CBANDS; b++ )
  1406.         for ( k = s3ind[b][0]; k <= s3ind[b][1]; k++ )
  1407.         dest[b] += s3_l[b][k] * source[k];    
  1408. }
  1409.